
/*
Output: Table 4. Observed and Predicted Percent Area Harvested
		Table SI1. Annual Percent Area Harvested, 1986-2019, Excluding Bedford County
		Table SI2. Annual Percent Area Harvested, 1986-2019, Including Bedford County
*/

run "/Users/holyfantastic/Dropbox/Forest/Paper/Nature/FinalDocuments/Code/000declare_path" // change to your local path

	
**# Housekeeping and clean data 
 

set scheme white_jet

use "ParcelLevelDataSept282023_VerNov82024_R1.dta" , clear

* demean ln_price
qui: sum ln_price
gen ln_price_dm = ln_price - r(mean)

* clean data
qui:{


la var forestAcres "Forested acres of parcel"
la var forestAcresTotal "Total forested acres owned by owner"

*Requiring at least one acre of forest
drop if forestAcres<1
drop if taxidnum == ""

*Focusing on pre enrollment years (pre-2020)
drop if year>2019

*Mean harvesting by year
bysort year: egen pcut_mean=mean(Landtr_Cut_mean_forest)


*clean variables
gen forestAcres2=forestAcres^2
gen distance2=distance^2
gen long_owner = YearsOwned > 10


*recode ndvi
replace ndvi_landtr_mean_forest = . if ndvi_landtr_mean_forest <=0 
replace ndvi_landtr_mean_forest = ndvi_landtr_mean_forest / 1000

*gen log of ndvi
gen ln_ndvi  = ln(ndvi_landtr_mean_forest * 1000 + 1)

*demean ndvi and ndvi2

foreach var in ndvi_landtr_mean_forest  {
	qui:sum `var'
	replace `var' = `var' - r(mean)
}
gen  ndvi_landtr_mean_forest2 = ndvi_landtr_mean_forest^2
gen  ndvi_landtr_mean_forest3 = ndvi_landtr_mean_forest^3


*Let's do everything with the percent harvested, not the share harvested
replace Landtr_Cut_mean_forest = Landtr_Cut_mean_forest * 100

**# regression analysis

**## look at key variables

tab group
egen taxidnum_id  = group(taxidnum)
cap egen county_id = group(county)
xtset  taxidnum_id year
sum Landtr_Cut_mean_forest
}

**## Table 4 - new
*c1 current owner
global parcel_controls  i.forestHaGroup i.roadGroup
global owner_controls  i.yearGroup i.iPerson ///
				i.iLocal
				
reghdfe Landtr_Cut_mean_forest  i.group  l.ndvi_landtr_mean_forest l.ndvi_landtr_mean_forest2  ///
											, a(i.county_id##i.year) vce(cl taxidnum_id year )

estadd local county_year_fe "Yes" , replace
estadd local parcel_controls "No" , replace
estadd local owner_specific "No" , replace
estadd local owner "All" , replace
qui: sum Landtr_Cut_mean_forest 
estadd local mean_y = round(r(mean),0.01)
estadd local method  = "One-step" , replace
estadd local N =  string(e(N), "%9.0f") , replace


eststo table_2_c1

				
				
*c2 Add parcel characteristics

reghdfe Landtr_Cut_mean_forest i.group l.ndvi_landtr_mean_forest l.ndvi_landtr_mean_forest2 $parcel_controls ///
											, a(i.county_id##i.year) vce(cl taxidnum_id year )

estadd local county_year_fe "Yes" , replace
estadd local parcel_controls "Yes" , replace
estadd local owner_specific "No" , replace
estadd local owner "All" , replace
qui: sum Landtr_Cut_mean_forest 
estadd local mean_y = round(r(mean),0.01)
estadd local method  = "One-step" , replace
estadd local N =  string(e(N), "%9.0f") , replace
eststo table_2_c2

*c3 Add owner characteristics and limit to current owners

reghdfe Landtr_Cut_mean_forest i.group l.ndvi_landtr_mean_forest l.ndvi_landtr_mean_forest2 ///
												$parcel_controls $owner_controls ///
												if year >= YearBought ///
											, a(i.county_id##i.year) vce(cl taxidnum_id year )

estadd local county_year_fe "Yes" , replace
estadd local parcel_controls "Yes" , replace
estadd local owner_specific "Yes" , replace
estadd local owner "Current" , replace
qui: sum Landtr_Cut_mean_forest if year >= YearBought
estadd local mean_y = round(r(mean),0.01)
estadd local method  = "One-step" , replace
estadd local N =  string(e(N), "%9.0f") , replace
eststo table_2_c3

				
foreach tab_path in table  {
	
	 	dis  "$`tab_path'/table_4_new.tex"

#delimit ; 
esttab table_2_c1 table_2_c2  table_2_c3
	using  "$`tab_path'/table_4_new.tex"
, replace f  compress
b(%12.3f) se(%12.3f) star(* 0.10 ** 0.05 *** 0.01)  
keep(1.group 2.group 3.group  L.ndvi_landtr_mean_forest L.ndvi_landtr_mean_forest2 )      coeflabels(L.ndvi_landtr_mean_forest " Lag.NDVI " L.ndvi_landtr_mean_forest2 " Lag.NDVI$^2$ ")
 label booktabs noobs nonotes  collabels(none) alignment(D{.}{.}{-1}) 
 stats( mean_y county_year_fe parcel_controls owner_specific owner  r2_a  N , labels( "Mean of Dep" "County-year FE" "Parcel controls" "Owner-specific controls" "Ownership"   "\hline Adjusted R-squared" "N" ))
nogaps nomtitles ;
 #delimit cr
 
	
}


**# Appendix: No Bedford
preserve

keep  if county != "Bedford"

*c1 current owner
global parcel_controls  i.forestHaGroup i.roadGroup
global owner_controls  i.yearGroup i.iPerson ///
				i.iLocal
				
reghdfe Landtr_Cut_mean_forest i.group l.ndvi_landtr_mean_forest l.ndvi_landtr_mean_forest2  ///
											, a(i.county_id##i.year) vce(cl taxidnum_id year )

estadd local county_year_fe "Yes" , replace
estadd local parcel_controls "No" , replace
estadd local owner_specific "No" , replace
estadd local owner "All" , replace
qui: sum Landtr_Cut_mean_forest 
estadd local mean_y = round(r(mean),0.01)
estadd local method  = "One-step" , replace
estadd local N =  string(e(N), "%9.0f") , replace
eststo table_2_c1

				
				
*c2 Add parcel characteristics

reghdfe Landtr_Cut_mean_forest i.group l.ndvi_landtr_mean_forest l.ndvi_landtr_mean_forest2 $parcel_controls ///
											, a(i.county_id##i.year) vce(cl taxidnum_id year )

estadd local county_year_fe "Yes" , replace
estadd local parcel_controls "Yes" , replace
estadd local owner_specific "No" , replace
estadd local owner "All" , replace
qui: sum Landtr_Cut_mean_forest 
estadd local mean_y = round(r(mean),0.01)
estadd local method  = "One-step" , replace
estadd local N =  string(e(N), "%9.0f") , replace
eststo table_2_c2

*c3 Add owner characteristics and limit to current owners

reghdfe Landtr_Cut_mean_forest i.group l.ndvi_landtr_mean_forest l.ndvi_landtr_mean_forest2 ///
												$parcel_controls $owner_controls ///
												if year >= YearBought ///
											, a(i.county_id##i.year) vce(cl taxidnum_id year )

estadd local county_year_fe "Yes" , replace
estadd local parcel_controls "Yes" , replace
estadd local owner_specific "Yes" , replace
estadd local owner "Current" , replace
qui: sum Landtr_Cut_mean_forest if year >= YearBought
estadd local mean_y = round(r(mean),0.01)
estadd local method  = "One-step" , replace
estadd local N =  string(e(N), "%9.0f") , replace
eststo table_2_c3

				
foreach tab_path in table  {
	
	 	dis  "$`tab_path'/table_4_new_app_R1.tex"

#delimit ; 
esttab table_2_c1 table_2_c2  table_2_c3
	using  "$`tab_path'/table_4_new_app_R1.tex"
, replace f  compress
b(%12.3f) se(%12.3f) star(* 0.10 ** 0.05 *** 0.01)  
keep(1.group 2.group 3.group  L.ndvi_landtr_mean_forest L.ndvi_landtr_mean_forest2 )      coeflabels(L.ndvi_landtr_mean_forest " NDVI " L.ndvi_landtr_mean_forest2 " NDVI$^2$ ")
 label booktabs noobs nonotes  collabels(none) alignment(D{.}{.}{-1}) 
 stats( mean_y county_year_fe parcel_controls owner_specific owner  r2_a  N , labels( "Mean of Dep" "County-year FE" "Parcel controls" "Owner-specific controls" "Ownership"   "\hline Adjusted R-squared" "N" ))
nogaps nomtitles ;
 #delimit cr
 
	
}

restore

**# Appendix: with Bedford 

reghdfe Landtr_Cut_mean_forest i.group l.ndvi_landtr_mean_forest l.ndvi_landtr_mean_forest2 ///
												$parcel_controls $owner_controls ///
												if year >= YearBought ///
												| county == "Bedford" ///
											, a(i.county_id##i.year) vce(cl taxidnum_id year )

			estadd local county_year_fe "Yes" , replace
			estadd local parcel_controls "Yes" , replace
			estadd local owner_specific "Yes" , replace
			estadd local owner "Current" , replace
			qui: sum Landtr_Cut_mean_forest if year >= YearBought
			estadd local mean_y = round(r(mean),0.01)
			estadd local method  = "One-step" , replace
			estadd local N =  string(e(N), "%9.0f") , replace
			
			eststo table_3_app_c1

reghdfe Landtr_Cut_mean_forest i.group l.ndvi_landtr_mean_forest l.ndvi_landtr_mean_forest2 ///
												$parcel_controls $owner_controls ///
												if year >= YearBought ///
												| (county == "Bedford" & (2019 - year < 20 )) ///
											, a(i.county_id##i.year) vce(cl taxidnum_id year )

			estadd local county_year_fe "Yes" , replace											
			estadd local parcel_controls "Yes" , replace
			estadd local owner_specific "Yes" , replace
			estadd local owner "Current" , replace
			qui: sum Landtr_Cut_mean_forest if year >= YearBought
			estadd local mean_y = round(r(mean),0.01)
			estadd local method  = "One-step" , replace
			estadd local N =  string(e(N), "%9.0f") , replace
			
			eststo table_3_app_c2
			
			
reghdfe Landtr_Cut_mean_forest i.group l.ndvi_landtr_mean_forest l.ndvi_landtr_mean_forest2 ///
												$parcel_controls $owner_controls ///
												if year >= YearBought ///
												| (county == "Bedford" & (2019 - year < 10 )) ///
											, a(i.county_id##i.year) vce(cl taxidnum_id year )

			estadd local county_year_fe "Yes" , replace											
			estadd local parcel_controls "Yes" , replace
			estadd local owner_specific "Yes" , replace
			estadd local owner "Current" , replace
			qui: sum Landtr_Cut_mean_forest if year >= YearBought
			estadd local mean_y = round(r(mean),0.01)
			estadd local method  = "One-step" , replace
			estadd local N =  string(e(N), "%9.0f") , replace
			
			eststo table_3_app_c3
			
			
foreach tab_path in table {
	
	 	dis  "$`tab_path'/table_4b_new_app_R1.tex"

#delimit ; 
esttab table_3_app_c1 table_3_app_c2  table_3_app_c3
, replace f  compress
b(%12.3f) se(%12.3f) star(* 0.10 ** 0.05 *** 0.01)  
keep(1.group 2.group 3.group  L.ndvi_landtr_mean_forest L.ndvi_landtr_mean_forest2 )      coeflabels(L.ndvi_landtr_mean_forest " Lag.NDVI " L.ndvi_landtr_mean_forest2 " Lag.NDVI$^2$ ")
 label booktabs noobs nonotes  collabels(none) alignment(D{.}{.}{-1}) 
 stats(  county_year_fe parcel_controls owner_specific owner  r2_a  N , labels(  "County-year FE" "Parcel controls" "Owner-specific controls" "Ownership"   "\hline Adjusted R-squared" "N" ))
nogaps 
mtitles("Add All Bedford" "Add Bedford, recent 20 years" "Add Bedford, recent 10 years")
 ;
 #delimit cr
 
	
}			

			
